Theory of Two-Dimensional Josephson Arrays in a Resonant Cavity 
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We consider the dynamics of a two-dimensional array of underdamped Josephson junctions placed 
in a single-mode resonant cavity. Starting from a well-defined model Hamiltonian, which includes 
the effects of driving current and dissipative coupling to a heat bath, we write down the Heisenberg 
equations of motion for the variables of the Josephson junction and the cavity mode, extending our 
previous one-dimensional model. In the limit of large numbers of photons, these equations can be 
expressed as coupled differential equations and can be solved numerically. The numerical results 
show many features similar to experiment. These include (i) self-induced resonant steps (SIRS's) 
at voltages V — nhfl/{2e), where Q is the cavity frequency, and n is generally an integer; (ii) a 
threshold number Nc of active rows of junctions above which the array is coherent; and (iii) a 
time-averaged cavity energy which is quadratic in the number of active junctions, when the array is 
above threshold. Some differences between the observed and calculated threshold behavior are also 
observed in the simulations and discussed. In two dimensions, we find a conspicuous polarization 
effect: if the cavity mode is polarized perpendicular to the direction of current injection in a square 
array, it does not couple to the array and there is no power radiated into the cavity. We speculate that 
the perpendicular polarization would couple to the array, in the presence of magnetic-field-induced 
frustration. Finally, when the array is biased on a SIRS, then, for given junction parameters, the 
power radiated into the array is found to vary as the square of the number of active junctions, 
consistent with expectations for a coherent radiation. For a given step, a two-dimensional array 
radiates much more energy into the cavity than does a one-dimensional array. 

PACS numbers: 05.45.Xt, 79.50.+r, 05.45.-a, 74.40.-(-k 



I. INTRODUCTION 

The properties of arrays of Josephson junctions have 
been of great interest for nearly twenty years Such ar- 
rays are excellent model systems in which to study such 
phenomena as phase transitions and quantum coherence 
in two dimensions. For example, if only the Josephson 
coupling energy is considered, the Hamiltonian of a two- 
dimensional (2D) array of Josephson junctions is formally 
identical to that of a 2D XY model [see, e. g. Ref. 
^. Arrays sometimes appear to mimic behavior seen in 
nominally homogeneous materials, such as high- Tc su- 
perconductors, which often behave as if they are com- 
posed of distinct superconducting regions linked together 
by Josephson coupling.^ Finally, the arrays are of poten- 
tially practical interest: they may be useful, for example, 
as sources of coherent microwave radiation if the indi- 
vidual junctions can be caused to oscillate in phase in a 
stable manner. 

Recently, our ability to achieve this kind of stable 
oscillation, and coherent microwave radiation, was sig- 
nificantly advanced by|-|a|rS|jjies of experiments by Bar- 
bara and collaborators .□'Ba'ElEl These workers placed two- 
dimensional underdamped Josephson arrays in a geome- 
try which allowed them to be coupled to a resonant mi- 
crowave cavity. The presence of the cavity caused the 
junctions to couple together far more efficiently than in 
its absence. As a result, the power radiated into the cav- 
ity has been found to be as much as 30% of the d. c. power 
injected into the array, far higher than the efficiency 
achieved in previous experiments. Even more surprising, 
this efhciency is achieved in underdamped arrays, which 



according to conventional wisdom should be especially 
difficult to synchronize, since each such junction exhibits 
bistability and hysteresis as a function of the external 
control parameters. These experiments haAt e . sJjrHul ated 
many theoretical attempts to explain them.clllj^EJE3 

In our previous work, we have presented a simple one- 
dimensional (ID) model which seems to account for npiW^ 
features of the observed cavity-induced coherence. tllEHl 
Despite the geometrical differences, the ID model does 
a surprisingly good job of capturing the physics of the 
experiments. However, a truly realistic test requires that 
the model be extended to a geometry closer to the ex- 
perimental one. In this paper, therefore, we present the 
necessary extension to 2D. Our results give significant in- 
sight into why the ID model works so well. In addition, 
they provide some clues about how one might understand 
experimental features which are still unexplained in ei- 
ther the ID or the 2D models. 

The remainder of this paper is organized as follows. 
In the next section, we describe our model Hamiltonian 
for a 2D current-driven, underdamped Josephson junc- 
tion array in a resonant cavity which supports a single 
mode. This Hamiltonian is a straightforward extension 
of that used in our previous work to describe ID arrays. 
In Section III, using this Hamiltonian, we write out the 
Heisenberg equations of motion for the junction variables 
and for the photon creation and annihilation operators 
for the cavity mode. We incorporate resistive dissipa- 
tion in the junctions in a standard way, by coupling the 
gauge-invariant phase differences across each junction to 
its own set of harmonic oscillator variables whose spec- 
tral density is chosen to produce Ohmic dissipation. In 
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the limit of large numbers of photons, we obtain classi- 
cal equations of motion for the variables. In Section IV, 
we present the numerical solutions of this model with an 
emphasis on features special to 2D, and we also give a 
comparison between the 2D and previous ID results. A 
concluding discussion and comparison with experiment 
follows in Section V. 
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II. MODEL HAMILTONIAN 

We will consider a 2D array of TV x M superconduct- 
ing grains placed in a resonant cavity, which we assume 
supports only a single photon mode of frequency f2. The 
array is thus made up of {N — 1)(M — 1) square plaque- 
ttes. There are a total of x Ny horizontal junctions, 
where = N — 1 and Ny = M. A current / is fed into 
each of the M grains on the left edge of the array, and 
extracted from each of the M grains on the right edge. 
Thus, the current is inject in the x direction, with no 
external current injected in the y direction. A sketch of 
this geometry is shown in Fig. |l|. We also introduce the 
terminology that a "row" of junctions, in this configura- 
tion, refers to a group of Ny junctions, all with left-hand 
end having the same x coordinate, and all being parallel 
to the bias current. One such row is indicated by the 
dashed lines in Fig. |l|. j— . 

In contrast to our previous work,E3 we will write the 
equations of motion for the grain variables (phases and 
charges) rather than junction variables, since in 2D, the 
junction variables cannot be treated as all independent 
(there are twice as many junctions as grains). 

We express our Hamiltonian in a form analogous to 
that of Ref. |l|: 

H = Hphoton + Hj + He + Hcurr + Hdiss- (1) 

Here Hphoton is the energy of the cavity mode, expressed 
as 

Hphoton = fi^ (^a^a + ^ , (2) 

where and a as the usual photon creation and anni- 
hilation operators. Hj is the Josephson coupling energy, 
and is assumed to take the form 

iI,, = -^£;,^-cos(7,,), (3) 

where Ef^ is the Josephson energy of the (ij)*'' junction, 
and 7,;j is the gauge- invariant phase difference across that 
junction (defined more precisely below). Efj is related to 
I^p the critical current of the (ij)*'* junction, by Ef^ — 
hlfj/q, where q = 2|e| is the Cooper pair charge. He is 
the capacitive energy of the array, which we write in a 
rather general form as 

Hc^lY.qHC-%n,n„ (4) 



FIG. 1: Sketch of the array geometry considered in our model. 
There are (Af x N) superconducting islands (black squares), 
making [{M — 1) x N + {N — 1) x M] Josephson junctions 
(crosses). An external current 7'^'^' is injected into each junc- 
tion at one end of the array and extracted from each junction 
at the other end. The array is placed in an electromagnetic 
cavity which supports a single resonant photon mode of fre- 
quency n. We have indicated by dashes a group of junctions 
which we denote a "row." Such a row is perpendicular to the 
current bias and is comprised of horizontal junctions. 

where C^^ is the inverse capacitance matrix, rii is the 
number of Cooper pairs on the z*'* grain, and g = 2e is 
the charge of a Cooper pair (we take e > 0). Note that in 
Ref. the variable rii was used to denote the difference 
between the numbers of Cooper pairs on the two grains 
comprising junction i. 

As in ID, the gauge-invariant phase difference, 7^, is 
the term which leads to coupling between the Josephson 
junctions and the cavity. We write it as 

- [(27r)/$o] / A-ds= cj), ~ cjjj - A,j , (5) 

J ij 

where (pi is the gauge-dependent phase of the supercon- 
ducting order parameter on grain i, and A is the vectar 
potential, which (in Gaussian units) takes the fornJijO 

A(x, t) = V{hc^)/m («W + «t W) E(x), (6) 

where E(x) is a vector proportional to the local electric 
field of the mode, normalized such that Jy (i^a;|E(x)p — 
1, and V is the cavity volume. The line integral is taken 
across the (ij)*'' junction. 

Given this representation for A, the phase factor Aij 
can be written 

Aij = gij{a -t- a^), (7) 
where gij takes the form 

I he? (27r)3 f , 

Clearly, gij is an effective coupling constant describing 
the interaction between the (ij)*'' junction and the cavity. 

We include a driving current and dissipation in a man- 
ner similar to that of Ref. ^ The driving current is 
included via a "washboard potential," Hcurr, of the form 

Hcurr — / T*7 5 (9) 

q ^-^ 

(*i>ll* 
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where I is the driving current injected in the x direc- 
tion into each grain on the left edge (and extracted 
from the right edge), and the sum runs over only those 
bonds in the x direction (each such bond is counted 
once). To introduce dissipation, each gauge- invariant 
phase difference, 7^, is coupled to a separate collec- 
tion ofpJiarpimiic oscillators with a suitable spectral 
density .liSlljO'Efl Thus, the dissipative term in the Hamil- 



tonian is 



H, 



diss 



diss 
ij 1 



(10) 



^0 = h/{Ae^) and Rij is a constant with dimensions of 
resistance (which proves to be the effective shunt resis- 
tance of the junction, as discussed below). 



III. EQUATIONS OF MOTION 



To obtain equations of motion, it is convenient to intro- 
duce the operators a = Ofl-t-ia/ and = aR — iai. These 
have the commutation relation [a^, a/] — i/2, which fol- 
lows [a, a^l = 1. In terms of these variables, 



where the sum runs over distinct bonds (ij), and 



fa,ij ^ij '^a,ij ~l" 



(j-^a.ij ) {_'^a,ij )^ 



{fa,i 



2 fT^a.ij ('^a,ij)^ 



{Pa,ij) 
^^^a.ii 



(11) 



The variables Ua^ij and Pa.ij, describing the a*'' os- 
cillator in the (ij)*'' junction, are canonically conjugate, 
and rua^ij and the mass and frequency of that 

oscillator. By choosing the spectral density, Jij{Lu), to be 
linear in \u!\, .wti assure that the dissipation in the junc- 
tion is ohmic.EjilZI We write such a linear spectral density 
as 



H photon — hQ.{a\ -\- of), 
and 7y takes the form 

lio = 0i - <t>j - 'igijCLR. 



(13) 



(14) 



The time-dependence of the various operators appear- 
ing in the Hamiltonian (|l|) is now obtained from the 
Heisenberg equations of motion. These are readily de- 
rived from the commutation relations for the various op- 
erators in the Hamiltonian (|l|). Besides the relations al- 
ready given, the only non-zero commutators are 



-ih SaJ3 Sij^k 



(15) 
(16) 



Jiji^) = IT- O^ij \^\ 0(^c -t^), 

Ztt 



(12) 



where Uc is a high-frequency cutoff (at which the as- 
sumption of ohmic dissipation begins to break down). 



e(wc 



is the usual step function, and aij is a dimen- 



sionless constant. We write it as atj = Ro/Rij, where 



where the last delta function vanishes unless (ij) and {k£) 
refer to the same junction. 

Using all these relations, we find, after a little algebra, 
the following equations of motion for the operators ^i, 
nj, aR, and a/: 



j 

jext 



d_R = a/, 
d/ = 



ft aR + '^ gij -j^ sin( 



Ef E(/. 



{fa,ij) 

a. 11 ^aAj ' n 



I a 



4>i ~ '^guaR) 



jext 

2gijai?) ^ g^j 




(17) 

(18) 
(19) 



(20) 



Here, the index I ranges over the nearest-neighbor grains the only external currents If^* are those along the left 
of i. In writing these equations, we have assumed that and right edges of the array, where they are i/'^^* [cf. 
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Fig. 1^. Eqs. (p7|)-(pO|) are equations of motion for the 
operators an, aj, rij, and (or 7^). In order to make 
these equations amenable to computation, we wiU later 
regacjd these operators as c- numbers, as we did earlier in 
1DJ13 This approximation is expected to he.reasonable 
when there are many photons in the cavity.c3 

The equations of motion for the harmonic oscillator 
variables can also be written out explicitly. However, 
since we have no direct interest in these variables, we 
instead eliminate them in order to incorporate a dissipa- 
tive term directly into the equations of motion for the 
other variables. Such a replacement is possible provided 
that the spectral density of each junction is linear in fre- 



quency, as noted above. In that case,ll3tall3'tZlt3 the os- 
cillator variables can be integrated out. The effect of 
carrying out this procedure is that one should make the 
replacement 



{fa,ij) 



% (21) 



wherever this sum appears in the equations of motion. 
Making the replacement (^) in Eqs. ( p^ ) and (pO|), and 
simplifying, we obtain the equations of motion for nj and 
ai with damping: 



ai 



iv) 



Text 



E l Rq ■ 
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jext 
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(«i>ll* 



(22) 
(23) 



Once, again, the index j is summed only over the nearest- 
neighbor grains of i. Equations (0), (|l|), (||), and (H) 
form a closed set of equations which can be solved for the 
time-dependent functions 7^, n^, and a/, given the ex- 
ternal current and the other parameters of the problem. 

It is now convenient to express these equations of mo- 
tion in terms of suitable scaled variables. We therefore in- 
troduce a dimensionless time r = tqRI'^/h = LOrt, where 
R and are suitable averages over Rij 
define the other scaled variables 



and Ifi . We also 



Rij — 

n = 
= 

aR,i = 



Rjj 
R ' 

n 

UJ 

I 

RI'=' 



o R 

2TT—aRj, 

Ko 



9ij = 

Cii — 



Rq 



27ri?" ■ 

UJ-j-RCij . 



(24) 
(25) 
(26) 
(27) 
(28) 

(29) 
(30) 



We assume that this takes the formi 

Cij = (Cd + Zi Cc) Sij 

— Cc (<^ij+x + <5ij'-x + ^i,j+y + ^i,j-y) ) (31) 

i. e., that there is a nonvanishing capacitance only 
between neighboring grains and between a grain and 
ground. Here Zi{— 4) is the number of nearest neigh- 
bors of grain i, Cd and Cc are respectively the diagonal 
(self) and nearest-neighbor capacitances, and x and y are 
unit vectors in the x and y directions. The correspond- 
ing Stewart-McCumber parameters are (3c — lOtRCc and 

Pd — LOrRCd- 

In Eq. (p7|), we introduced the potential Vi on site i, 
which is expressed through the number variables rij as 



lY^iC ^)ijnj. 



(32) 



The integral of the electric field across junction (ij) is 
written in terms of the Vi 's as 



Vi - V-j - 2gijilaj. 



(33) 



The last equation involves the capacitance matrix Ci- 



Carrying out these variable changes, we find, after 
some algebra, that the equations of motion can be ex- 
pressed in the following dimensionless form: 



(34) 
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ar 



J L / \ 



3 



sin(0j - (pi - 2gjiaR) + -^{Vi - Vi - 2guClai) 



R 



iij sm{4>i-(l)j-2gijdR) + -^{Vi - Vj - 2Clgijdi) 



{u>ll* 



(35) 

(36) 
(37) 



These equations are readily generalized to treat external 
currents with non-zero components in both the x and the 
y directions, and to geometries other than lattices with 
square primitive cells. 



[cf., e. g., Refs. ^ and In our simulations for 2D 
arrays coupled to a resonant cavity, we observe this same 
phenomenon, as discussed below. 



IV. NUMERICAL RESULTS 

We solve Eqs. (|3^) ~ (|3^) numerically,r-by implement- 
ing the adaptive Bulrisch-Stoer method,^ as described 
further in Ref. For simplicity, we assume that the 
coupling constants gij have only two possible values, gx 
and corresponding to junctions in the x and y direc- 
tion respectivelyB3 This assumption should be reasonable 
if two conditions are satisfied: (i) there is not much disor- 
der in the characteristics of the individual junctions; and 
(ii) the wavelength of the resonant mode is large com- 
pared to the array dimensions. Although assumption (ii) 
is not obviously satisfied for the experimental arrays, the 
model may still be reasonable in certain array/cavity ge- 
ometries, as discussed further below. 

In the presence of a vector potential, it is customary 
to definje-a frustration (t) for the plaquette by the 
relatioHCj 

plaquette plaquette 

where the sum runs over bonds in the /i*'* plaquette. For 
a general position-dependent cnj, /^(t) ^ 0, but if gx and 
(jy are both position-independent, then ff_i{T) = 0. The 
possibility of having distinct coupling constants gx and 
gy along x and y bonds arises from differences in possible 
polarizations of the resonant mode, and leads to effects 
which cannot be captured in a ID model, as discussed 
below. 

Before discussing our numerical results, we briefly sum- 
marize one well-known feature of underdamped Joseph- 
son arrays in the absence of coupling to a resonant cavity. 
At certain applied currents, the individual junctions in 
such an array are bistable - that is, they can be placed 
in an "active" (resistive) or an "inactive" (superconduct- 
ing) state, by a careful choice of initial conditions. For 
an applied current in the x direction, when a single hor- 
izontal junction is chosen to be in the active state, it is 
found that all the other horizontal junctions in the same 
"row" (cf. Fig. ^ also go active, provided that there is 
at least a little disorder in the junction critical currents 



A. Horizontal coupling 

We first consider the case gx 9y — 0, with driving 
current parallel to the x axis. In Fig. ^, we show a se- 
ries of current-voltage (IV) characteristics for this case. 
We consider an array of 10 x 4 grains, with capacitances 
/?c = 20 and = 0.05, gx = 0.012, and Q = 0.41. The 
critical current through the (ij)*'' junction is /f^ = 1-1- 
where the disorder Ay is randomly selected with uni- 
form probability from [—A, A]. In this plot, A = 0.05. 
The product IfjRij is assumed to be the same for all 
junctions, in, accordance with the Ambegaokar-Baratoff 
expression.Ej In addition, (3^ and Pc are assumed to be 
the same for all junctions. The calculated IV's are shown 
as a series of points. The directions of the arrows indicate 
whether the curves were obtained under increasing or de- 
creasing current drive, or both. The horizontal dashed 
curves correspond to voltages where self-induced resonant 
steps (SIRS's) are expected, namely {V)r/{NRIc) = 
in our units, where {V)r denotes the time-averaged volt- 
age. The dotted lines are guides to the eye. Each nearly 
horizontal series of points denotes a calculated IV char- 
acteristic for a different number of active rows Na, and 
represents Na x iV^ (horizontal) junctions sitting on the 
first integer (n = 1) SIRS. The calculated voltages for the 
various Na 's agree well with the expected values given by 
the dashed horizontal lines. The long straight diagonal 
line segment, which is common to all the different A'^'s, 
represents the ohmic part of the IV characteristic with 
all rows active. For the sake of clarity, we have chosen 
not to plot the corresponding segments for other choices 
of Na < 10. Besides the integer SIRS's we find that for 
this 2D array, it is possible to bias individual active rows 
on either the n = 1/2 or the n = 2 SIRS. (A small seg- 
ment of an n = 1/2 case is visible in the lower left of the 
figure.) Similar behavior is found in the case of Shapiro 
steps produced by a combined d. c. and a. c. current in 
a conventional underdamped Josephson junction (see, e. 
g. Ref. 13). 

Although the full hysteresis loop is shown in Fig. |^ only 
for Na — 10 active rows, the IV curves for other values 
of Na are also hysteretic. Specifically (as also found pre- 
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FIG. 2: Calculated current-voltage characteristic for a 10 x 4 
array with cavity frequency Q = 0.41, capacitance parame- 
ters /3c = 20 and f3d = 0.05, disorder parameter A = 0.05 
and junction-cavity coupling in the horizontal direction Qx ~ 
0.012. The horizontal dashed lines show voltages at which 
the various SIRS's are expected. These correspond to differ- 
ent numbers of rows of horizontal junctions in the active state. 
Arrows denote that the given IV was taken in the direction 
of increasing or decreasing current. 



viously in the ID case), whenever Na > 4, the number of 
active rows increases when the SIRS's become unstable. 
That is, if the current is increased so that a given SIRS 
becomes unstable, the IV characteristic jumps up onto a 
higher SIRS, and also some of the individual rows jump 
onto the n — 2 SIRS. The IV curve only jumps onto the 
ohmic branch if ///c > 1. By contrast, if the applied cur- 
rent is changed so that the SIRS's become unstable for 
Na ^ 4, the number of active rows remains unchanged 
and the IV curve immediately becomes ohmic. In this 
regime, if / is increased so that I/Ic ~ 1, all the re- 
maining horizontal junctions become active and the IV 
characteristic also becomes ohmic. Another feature of 
these results worth noticing is that the width of the SIRS 
plateaus is non-monotonic in Na- By "width" of an SIRS, 
we mean the range of driving currents for which the SIRS 
is stable. 

Fig. H shows the IV characteristics for three differ- 
ent arrays, each with all rows in the active state: (i) 
a 40 X 1 (full curve), (ii) a 40 x 2 (dotted curve) and 
(iii) a 40 X 3 (long-dashed curve). Each array has the 
parameters = 0.015, Cl = 0.49, (3c = 20, /3d = 0.05 
and A — 0.05. Once again, the arrows denote the direc- 
tions of current sweep. The horizontal dot-dashed curve 
shows the expected position of the SIRS corresponding 
to Na = 40 [V/{N^RIc) = fi]. The curves show that all 
three arrays have qualitatively similar behavior. First, if 
the array is started from a random initial phase configu- 
ration, such that / = I/Ic > 1 -I- A, and / is decreased, 



FIG. 3: Calculated current-voltage characteristics for a 40 x 1 
(full hue) , a 40 x 2 (dotted line) and a 40 x 3 (long-dashed line) 
array, all with parameters = 0.015, Q. = 0.49, /3c = 20, 
/3d — 0.05 and A = 0.05. The horizontal dot-dashed line 
shows the expected position of the SIRS. Note that as the 
array width increases, the smallest value of / at which all the 
active junctions phase-lock on the SIRS also increases, and 
IV characteristic on the SIRS has an increasing bend. Hence, 
increasing the array width at fixed Qx has a effect similar to 
that of increasing Qx at fixed width. The arrows indicate the 
direction of the current sweep. 



then all the rows lock on to the Na = 40 SIRS. Secondly, 
if / is further decreased, the Na = 40 active state eventu- 
ally becomes unstable and all the junctions go into their 
superconducting states. Finally, if / is increased start- 
ing from a state in which the array is on the Na = 40 
SIRS, the SIRS remains stable until / reaches the critical 
current for the various rows, and the IV curve becomes 
ohmic. 

The behavior shown in Fig. |^ with increasing array 
width is very similar to that found previously in ID ar- 
rays with increasing coupling strength. In other words, 
the key parameter in understanding the curves of Fig. 

is the product Nyg^- For example. Fig. || shows that 
the effect of increasing Ny while keeping g^ constant is 
to raise slightly the maximum value of / for which the 
active junctions are still locked onto the Na = 40 SIRS. 
Furthermore, that portion of the SIRS which corresponds 
to small / is not perfectly flat (i. e. not at the expected 
constant voltage V / {N^RIc) — but instead increases 
slightly with increasing / [cf. Fig. ||. The degree of this 
non-flatness increases with increasing Ny. Precisely anal- 
ogous effects jaxe seen in calculations for ID arrays with 
increasing cjxt^ This is another piece of evidence that the 
key parameter is the product NyCjx. 

In Fig. ^ we plot the time-averaged energy E{Na) = 
(a|j + Cj),- in the cavity for three different arrays: 40 x 1 
(stars), 40 X 2 (circles), and 40 x 3 (squares). In all cases. 
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FIG. 4: Time- averaged scaled energy E in the resonant cavity 
as function of active number of rows for a 40 x 1 (asterisks), 
a 40 X 2 (circles) and a 40 x 3 (squares) array with driving 
current I — 0.58. All the other parameters are the same as 
those of Fig. ^. Inset: an enlargement of the IV characteristics 
near the synchronization threshold, on a logarithmic vertical 
scale. Note that the threshold number of active junctions for 
synchronization decreases with increasing array width. 



/ = 0.58, and the other parameters are the same as those 
of Fig. ||. Below a threshold value of A^o, (which we de- 
note Nc and which depends on A^j,), the active rows are 
in the McCumber state (not on the SIRS's). In this case, 
E{Na) is small and shows no obvious functional depen- 
dence on Na (see inset). By contrast, above threshold, 
E{Na) is much larger and increases as N^. 

Fig. ^ shows that, when Ny is increased at fixed cjx, Nc 
decreases. Precisely this same trend is observed when we 
increase gx while holding Ny fixed (and was observed in 
our previous ID calculations with increasing g^). Thus, 
once again, the relevant parameter in understanding the 
threshold behavior appears to be Nycj^. 

As in ID arrays, it is useful to introduce a Kuramoto 
order parameter which describes the phase ordering. For 
the 2D arrays, we define a Kuramoto order parameter 
(r/,,)^. for the horizontal bonds by 



1 



A^aA^^ 



(I E ^'"'^1 



(39) 



where Na is the number of active rows, A'j^ is the number 
of horizontal junctions in a single row and the sum runs 
over all the active, horizontal junctions. [The analogous 
quantity (r„),- for the vertical junctions is irrelevant when 
gv = 0, since in this case these junctions are inactive.] For 
the parameters shown in Fig. ^, we have found, as in our 
previous ID calculations, that {rh)T ^ 1 for Na > Nc 
while {rh)T <C 1 for Na < Nc- This behavior (which we 
do not show in a figure) reflects the fact that, for the 
value of / used in Fig. 0, none of the active junctions are 



FIG. 5: The time-averaged Kuramoto order parameter (rh)^ 
(defined in Eq. ( |39| ) as a function of the number of active 
rows on SIRS's for a 20 x 2 Josephson array with Cl — 0.49, 
gj: = O.OI, /3c = 20., f3d = 0.05, A = 0.1 and bias current 
/ = 0.53. 



on a SIRS when Na < Nc; hence, these junctions are not 
in phase with one another, and the value of {rh)T reflects 
this lack of coherence. 

For certain array parameters, / can be chosen so that 
all the active junctions lie on SIRS's, however many ac- 
tive rows Na there are. [In Fig. ||, for example, J ~ 0.5 
would achieve this result.] In such cases, even though 
all the active junctions are oscillating with the same fre- 
quency, and locked onto SIRS's, it is still possible to have 
(?'h)T < 1- In this situation, the Kuramoto order param- 
eter {rh,n)T ~ 1 for the individual rows. This occurs 
because the rows are not perfectly phase-locked to one 
other. An example of such behavior is shown in Fig. ||, 
for a 20 X 2 junction array for several numbers Na of ac- 
tive rows. The other parameters are Ct = 0.49, gx = 0.01, 
/3c = 20., l3d = 0.05, A = 0.1 and / = 0.53. As the num- 
ber of active rows on the SIRS's increases, {rx)T 1- 
[Also, of course, {rx)T = 1 for one active row on a SIRS.] 
Numerically, we find that it is easier in 2D than in ID 
to achieve a state with all active junctions biased on a 
SIRS, but with {rx)T < 1- In all such cases, we can easily 
cause {rx)T ^ 1 simply by increasing gx- 

The threshold shown in Fig. ^ corresponds to a tran- 
sition from a state in which none of the active junctions 
are on SIRS's to a state in which all are on SIRS's. It 
is possible to choose / so as to have any number of ac- 
tive rows Na on SIRS's. In this case, the cavity energy 
E(Na), in our model, is approximately quadratic in Na, 
with no obvious threshold behavior. This feature of our 
results is discussed further below. 
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FIG. 6; Phase plot of points (7, 7) where 7 is the gauge-invariant phase difference across a junction, shown for (a) a vertical 
junction and (b) a horizontal junction, in an array in which the cavity couples only to the horizontal junctions: = 0., 
Qy = 0.5. The array size is 10 x 4, and the other parameters are Q, = 0.45, A = 0.05, (3c = 20, f3d = 0.05 and / = 0.46. The 
vertical junction in (a) displays aperiodic motion with very small amplitude, corresponding to no time-averaged voltage drop 
across that junction, while the horizontal junction in (b) has a phase difference which varies periodically in time. 



B. Vertical coupling 

We have also investigated the case of gx — 0, (jy ^ 0, 
for a wide range range of gy values. For our geometry, we 
have not been able to find any value for gy for which a 
SIRS develops. In essence, when the cavity couples only 
to the vertical junctions, it is invisible in the IV char- 
acteristics. This behavior is easily understood. In this 
geometry, with current applied in the x direction, both 
the time-averaged voltage and the time-averaged current 
through the vertical junctions are very small. Hence, too 
little power is dissipated in the vertical junctions to in- 
duce a resonance with the cavity. 

To illustrate this behavior, we show in Fig. ^ some 
representative phase plots of ("fij , jij ) for (a) a vertical 
junction and (b) a horizontal junction in a 10 x 4 array 
with gx ^ 0, gy ^ 0.5, (3c = 20, Pa = 0.05, Q = 0.45, 
A = 0.05 at bias current / = 0.46 (close to a possible 
resonance with cavity). The phase plot for the verti- 
cal junction exhibits small-amplitude aperiodic motion, 
while that of the horizontal junction shows that this junc- 
tion is in its active state and undergoing periodic motion 
in phase space. This lack of response by the y junctions 
to the cavity probably explains why the ID simulations 
describe the experiments so well. 

It is no surprise that the cavity interacts only very 
weakly with the vertical junctions. From previous stud- 
ies of both underdamped and overdamped disordered 
Josephson arrays in a rectangular geometry (see, e. g., 
Refs. and[25|), it is known that when current is applied 
in the x direction, the y junctions remain superconduct- 
ing, with {V)r ~ 0, while the x junctions comprising 
an active row are almost perfectly synchronized, with 



(r.) « 1. 

If there were an an external magnetic field perpendic- 
ular to the array, we believe that SIRS's would be gen- 
erated for gy 7^ 0, even if g^ = 0. In this case, the mag- 
netic field would induce frustration. Specifically, since 
the sum of the gauge invariant phase differences around a 
plaquette must be an integer multiple of 2tt, the presence 
of magnetic-field-induced vortices piercing the plaquettes 
would induce nonzero voltages across, and supercurrents 
in, the y junctions. It would be of great interest if calcu- 
lations were carried out in such applied magnetic fields. 



C. Comparison with ID Model 

We now compare our 2D results explicitly with those 
for ID arrai/s. In our earlier ID model, we found 
numericallytj that the threshold number of active junc- 
tions Nc was inversely proportional to the coupling con- 
stant g. This behavior is reasonable because the inho- 
mogeneous term driving the cavity variable cir is pro- 
portional to the product of g and Na- 

Some of our numerical trends in the 2D case can be 
understood similarly. For example, the inhomogeneous 
term in Eq. ( p7| ) is the last term on the right-hand side. 
It is proportional to the sum of the coupling constants 
gij over all the junctions parallel to I. Thus, for g^ ^ 0, 
(jy = 0, and for the same driving current /, we expect 
that an x Ny array with a coupling constant g^ should 
behave like an x 1 array with coupling constant Nyg^. 

To check this hypothesis, we compare, in Fig. |^, the 
IV characteristics of a 10 x 1 array having coupling con- 
stant gx;Wy.i = 0.0259 with those of a 10 x 10 array with 
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Current I/I, Current I/I, 



FIG. 7: IV characteristics for a 10 x 1 array (*) and a 10 x 10 
array (o). The 10 x 1) array has parameters gx,ioxi = 0.0259, 
D. = 0.41, /3c = 20, = 0.05 and A = 0.05. The expected 
position of the SIRS's are marked by horizontal dashed lines. 
The 10 x 10 array has gi:,ioxio = 0.00259, and the other pa- 
rameters are the same as for the 10 x 1 array. The IV charac- 
teristics are shown for both increasing and decreasing current 
drive, as discussed in the text. 



coupling constant gx-.ioxw = 0.00259. The other param- 
eters are the same for the two arrays: = 0.41, f3c = 20, 
Pd = 0.05 and A — 0.05. The expected positions of the 
SIRS's [at V/{NRIc) = are indicated by dashed hor- 
izontal lines. Indeed, the two sets of IV characteristics 
are very similar. Even some of the subtle differences can 
be understood in a simple way. For example, the 10 x 10 
IV's are slightly flatter than the 10x1 curves. We believe 
this extra flatness occurs because the individual junction 
couplings in the 10 x 1 array are 10 times larger than 
those in the 10 x 10 array. From our previous ID sim- 
ulations, the IV's on the steps become more and more 
rounded as gx increases, i. e., the voltage on the lower 
portion of the SIRS is no longer independent of / [cf. Ref. 
|l2| . Precisely this behavior is seen in Fig. ^ 

Another subtle difference between the ID and 2D 
curves of Fig. |^ is the values of the so-called "retrap- 
ping current" in the two sets of curves (i. e. the current 
values below which the McCumber curve becomes unsta- 
ble). We beheve that this difference can be understood 
in terms of the effects of disorder in the junction critical 
currents in ID and 2D. Speciflcally, for a given value of 
A, the 2D arrays are effectively less disordered than the 
ID arrays, since the average critical current for a single 
row has a smaller rms spread than the critical current of 
a single junction in a ID array. 

An important similarity between the two sets of curves 
is that, in both the 10 x 10 and the 10x1 arrays, the width 
of the SIRS's varies similarly (and non-monotonically) 
with the number of active rows. This behavior distin- 



FIG. 8: Time-averaged reduced cavity energy E, for a 10 x 
1 array and a 10 x 10 array for the same choice of array 
parameters as in Fig. ^. The calculations are carried out on 
the decreasing current branch with all rows active. Note that 
Qx for the 10 x 10 array is 10 times smaller than that of the 
10 X 1 array. 



guishes our predictions from some other models ,aE3 in 
which the cavity is modeled as an RLC oscillator con- 
nected in parallel to the entire array, and whichjpredicts 
a monotonic dependence of SIRS width on /Vq.eI 

In Fig. H we plot the reduced time-averaged cavity en- 
ergy E — {a\ + a'j)r as a function of J = I/Ic for both 
arrays of Fig. |^, under conditions such that all rows are 
active. This plot is obtained by following the decreas- 
ing current branch. Surprisingly, when the 10 x 10 array 
(with .9(10x10) = 0.1 5(10x1)) locks on to the SIRS, E 
jumps to a value which is approximately two orders of 
magnitude larger than that of the corresponding jump 
in the 10 x 1 array, even though the parameter Nyg^ is 
the same for both arrays. We believe that the difference 
is due simply to the greater number of junctions which 
are driving the cavity in the 2D case. Even though the 
width of the steps is controlled primarily by the parame- 
ter NyCix, the energy in the cavity is determined by the 
square of the number of radiating junctions. This square 
is 100 times larger for the 2D array than for the ID array. 



V. SUMMARY 

In this paper, we have derived equations of motions 
for a 2D array of underdamped Josephson junctions in 
a single-mode resonant cavity, starting from a suitable 
model Hamiltonian and including the effects of both a 
current drive and resistive dissipation. In the limit of 
zero junction-cavity coupling, these equations of motion 
correctly reduce to those describing a 2D array of resis- 
tively and capacitively shunted Josephson junctions. 
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As in our previous ID model, the present equations 
of motion lead to a transition from incoherence to co- 
herence, as a function of the number of active rows Na- 
This transition again results from the effectively mean- 
ficld-likc nature of the interaction between the junctions 
and the cavity. Specifically, because each junction is, in 
effect, coupled to every other active junction via the cav- 
ity, the strength of the effective coupling is proportional 
to the number of active junctions. Thus, for any g^, no 
matter how small, a transition to coherence is to be ex- 
pected for sufficiently large number of active rows Na- 
We also found a striking effect of polarization: the tran- 
sition to coherence occurs only when the cavity mode is 
polarized so that its electric field has a component par- 
allel to the direction of current flow. 

Next, we briefly compare our timnerical results to 
the behavior seen in experiments.LHj Our calculations 
show the following features seen in experiments: (i) self- 
induced resonant steps (SIRS's) in the IV characteris- 
tics; (ii) a transition from incoherence to coherence above 
a threshold number of active junctions; and (iii) a to- 
tal energy in the cavity which varies quadratically with 
the number of active junctions when those junctions are 
locked onto SIRS's. There may, however, be some differ- 
ences as well. In particular, our transition to a quadratic 
behavior occurs when the active junctions are locked onto 
SIRS's. In possible contrast to our results, in some ex- 
perimental arrays,Ll it has been reported that even be- 
low the "coherence threshold," individual rows of junc- 
tions are locked onto SIRS's, but these SIRS's are not 
coherent with one another, and hence, do not radiate 
an amount of power into the cavity proportional to the 
square of the number of junctions on the SIRS's. Thus 
far, in our calculations, we have found that when Na junc- 
tions are locked onto the steps, the energy in the cavity 
is quadratic in Na- The threshold, in our calculations, 
occurs when all the active junctions lock onto SIRS's, 
not when active junctions which are already locked onto 
SIRS's become coherent with one another. 

For some choices of the parameters gx, ^, A, /3, and 
/, we find dynamical states such that all active rows lock 
onto SIRS's while {r)r < 1. In such states, the Kuramoto 
order parameter for the individual rows is still {r)r ~ 1, 
implying that the rows are not perfectly phase-locked 
to each other. An example of such a state is shown in 
Fig. ||. In such states, our calculated energy E in the 
cavity appears to vary smoothly with Na, and exhibits no 



threshold behavior, in contrast to what we find at other 
applied currents [cf. Fig. This behavior appears to 
differ from what was reported experimentally in a recent 
paper;! the reasons for the difference are not clear to us. 

In summary, we have extended our previous theory of 
Josephson junction arrays coupled to a resonant cavity 
to the case of 2D arrays. The 2D theory bears many 
similarities to the ID case, and makes clear why the ID 
model works so well. These similarities arise because, 
in a square array, the coupling to the cavity takes place 
only through those junctions which are parallel to the ap- 
plied current. Again as in ID, our model leads to clearly 
defined SIRS's whose voltages are proportional to the res- 
onant frequency of the cavity. Another similarity is that, 
in 2D as in ID, when a fixed number of rows are biased 
on a SIRS, the cavity energy is linear in the input power. 

We also find some results which are specific to 2D. 
For example, whenever one junction in a given row is 
biased on a SIRS, all the junctions in that row phase- 
lock onto that same SIRS. In addition, the time-averaged 
energy contained in the resonant cavity is quadratic in 
the number of active rows, but, when the array is biased 
on a SIRS, is much larger in the 2D array than in the 
ID array, for the same value of the coupling parameter 
9xNy. 

When the cavity mode is polarized perpendicular to 
the direction of current drive, we find that the cavity 
does not affect the array IV characteristics. Our equa- 
tions suggest that this non-effect might change if the ar- 
ray were frustrated, e. g., by an external magnetic field 
normal to the array. Such frustration would cause junc- 
tions in the x and y direction to be coupled. It would 
be of great interest if this speculation could be tested 
experimentally. 



Acknowledgments 

This work has been supported by the National Science 
Foundation, through grant DMROl-04987, and in part 
by the U.S. -Israel Binational Science Foundation. Some 
of the calculations were carried out using the facilities of 
the Ohio Supercomputer Center, with the help of a grant 
of time. We are very grateful for valuable conversations 
with Profs. T. R. Lemberger and C. J. Lobb. 



Electronic address: Almaas.l@nd.edu; Present address 



Deptartment of Physics, University of Notre Dame, Notre 
Dame, Indiana 46556, 



Electronic address: 3troud@mps.ohio-state.edu 



For an extensive review, see, e. g. R. S. Newrock, C. 
J. Lobb, U. Geigenmiiller, and M. Octavio, Solid State 
Physics 54, 263 (2000). 

See, for example, P. M. Chaikin and T. C. Lubensky. 



Principles of Condensed Matter Physics. Cambridge Univ. 
Press, New York, 1995. 

^ M. Tinkham. Introduction to Superconductivity. McGraw- 
Hill, New York, Second edition, 1996. 

* P. Barbara, A. B. Cawthorne, S. V. Shitov, and C. J. Lobb. 
Phys. Rev. Lett. 82, 1963 (1999). 

5 P. Barbara, G. Filatrella, C. J. Lobb, and N. F. Peder- 
sen, "Cavity Synchronization of Underdamped Josephson- 



11 



Junction Arrays," to bo published in Studies in High- 
Temperature Superconductors, edited by V. Narlikar 
(NOVA Science Publishers, New York, in press). 
B. Vasilic, P. Barbara, S. V. Shitov, and C. J. Lobb. 
Published in the proceedings of Applied Superconductivity 
Conference 2000. 
^ B. Vasilic, P. Barbara, S. V. Shitov, and C. J. Lobb. Phys. 

Rev. B65, 180503(R) (2002). 
* B. Vasilic, P. Barbara, S. V. Shitov, E. Ott, T. M. Anton- 
sen, and C. .]. Lobb. Abstract Y27.001 of the APS March 
Meeting, Seattle, WA., March 2001. 
® G. Filatrella, N. F. Pedersen, and K. Wiesenfeld. Phys. 
Rev. E61, 2513 (2000). 

1° G. Filatrella, N. F. Pedersen, and K. Wiesenfeld. IEEE 
Trans. Appl. Supercond. 

" E. Almaas and D. Stroud. Phys. Rev. B63, 144522 (2001). 
See also E. Almaas and D. Stroud, Phys. Rev. B 64, 
179902(E) (2001). 

1^ E. Almaas and D. Stroud. Phys. Rev. B65, 134502 (2002). 
J. C. Slater. Microwave Electronics. D. Van Nostrand, 
New York, 1950. 

A. Yariv. Quantum Electronics. J. Wiley & Sons, New 
York, Second edition, 1975. 

V. Ambegaokar, U. Eckern, and G. Schon. Phys. Rev. Lett. 
48, 1745 (1982). 



A. O. Caldeira and A. J. Leggett. Phys. Rev. Lett. 46, 211 
(1981). 

" A. O. Caldeira and A. J. Leggett. Ann. Phys. (N. Y.) 149, 
374 (1983). 

S. Chalcravaxty, G.-L. Ingold, S. Kivelson, and A. Luther. 

Phys. Rev. Lett. 56, 2303 (1986). 
1^ R. Fazio and G. Schon. Phys. Rev. B 43, 5307 (1990). 
^° B. J. Kim and M. Y. Choi. Phys. Rev. B 52, 3624 (1995). 
2^ W. H. Press, S. A. Teukolsky, W. T. Vetterling, and B. P. 

Flanncry. Numerical Recipes. Cambridge Univ. Press, New 



York, 1992. 

Note that one always has the relation gij = —Qji. 



22 

See, for example, S. Teitel and C. Jayaprakash, Phys. Rev. 
B27, 598 (1983). 

T. S. Tighe, A. T. Johnson, and M. Tinkham. Phys. Rev. 
B44, 10286 (1991). 
2^ W. Yu and D. Stroud. Phys. Rev. B 46, 14005 (1992). 
V. Ambegaokar and A. Baratoff, Phys. Rev. Lett. 10, 486 
(1963). 

J. R. Waldram. Superconductivity of Metals and Cuprates. 
Inst, of Phys. Publ., Philadelphia, 1996. 
^® W. Yu. Topics in the Dynamical Properties of Josephson 
Networks. PhD thesis, The Ohio State University, 1994. 



